%--------------------------------------------------------------------------
% computes the neutral stability curves 
%--------------------------------------------------------------------------


function [kk, M, fun] = lo_neutral_stab_curves(kk, figs)

fun = @(k) 0.24e2 * k ^ 3 * sinh(k) * (-cosh(k) * sinh(k) + k) / ((-0.3e1 * k - 0.6e1 * k ^ 3) * cosh(k) ^ 3 + 0.3e1 * cosh(k) ^ 2 * sinh(k) + (0.6e1 * k ^ 3 * sinh(k) ^ 2 + 0.3e1 * k + 0.3e1 * k ^ 3) * cosh(k) + 0.2e1 * (k ^ 4 - 0.3e1 / 0.2e1 + 0.3e1 / 0.2e1 * k ^ 2) * sinh(k));

for j = 1:length(kk)
    
    k = kk(j);
    
    if (k < 1e-12)
        M(j) = 80 + 128 / 63 * k^2;
    else
        M(j) = fun(k);
    end
    
end


if (figs)
    plot(kk, M, 'k','linewidth',1);
    xlabel('$\check{k}$','interpreter','latex','fontsize',12);
    ylabel('$\mathcal{M}$', 'interpreter','latex','fontsize',12);
end